#include "linbox/linbox-config.h"
#include "linbox/util/timer.h"
#include "linbox/field/modular.h"
#include "linbox/matrix/blas-matrix.h"
#include "linbox/algorithms/matrix-hom.h"
#include "gmp++/gmp++.h"

#define _LB_ITERS 6

using namespace LinBox ;

long modit (Integer & n, const long & l)
{

	// if (!n) return 0UL; /* trop cher */
	// bool isneg = (n<0) ;
	// CONDITION: mpz_tdiv_ui does NOT consider the sign of gmp_rep
	long res = mpz_tdiv_ui( (mpz_ptr)n.get_mpz(), l);
	if (!res) return 0UL ;
	if (n<0)
		res = l-res ;
	return res ;

}

	template<class T>
void clear (BlasMatrix<T> & A)
{

	typename BlasMatrix<T>::RawIterator it = A.rawBegin();
	for ( ; it != A.rawEnd(); ++it )
		*it = 0 ;
	return ;
}

	template<class T>
bool equal (const BlasMatrix<T> & A, const BlasMatrix<T> & B)
{
	if (A.getPointer() == B.getPointer()){
		std::cout << "bleuf !" << std::endl;
		return true;
	}

	typename BlasMatrix<T>::ConstRawIterator it = A.rawBegin();
	typename BlasMatrix<T>::ConstRawIterator jt = B.rawBegin();
	for ( ; it != A.rawEnd(); ++it,++jt ){
		// std::cout << *it << ',' << *jt << std::endl;
		if (*it != *jt) return false ;
	}
	return true ;
}

int main()
{
	int n    = 2000;   //taille des matrices
	int l    = 70;     //longueur des entiers
	double p = 65531 ; // caratéristique
	typedef Modular<double> Field;
	typedef typename Field::Element Element;
	Field F(p);
	BlasMatrix<Element>  Ap(n,n);
	BlasMatrix<Integer>  A(n,n);
	{ // random init
		typename BlasMatrix<Integer>::RawIterator it = A.rawBegin();

		for ( ; it != A.rawEnd(); ++it ) {
			*it = Integer::random<false>(l);
		}
	}

	typename BlasMatrix<Integer>::RawIterator it = A.rawBegin();
	Timer tim ;
	///////////////////////////////
	tim.clear();
	tim.start();
	for (size_t i = _LB_ITERS ; --i ; )  // répétitions
		MatrixHom::map(Ap,A,F);
	tim.stop();
	// A.write(std::cout << "--------A---------" << std::endl);
	// Ap.write(std::cout << "--------Ap--------" << std::endl);
	std::cout << "Hom   : " << tim << std::endl;
	const BlasMatrix<Element> Bp = (const BlasMatrix<Element>&) Ap ; // référence (oupa)
	///////////////////////////////
	clear(Ap);
	tim.clear();
	tim.start();
	typename BlasMatrix<Element>::RawIterator jt ;
	for (size_t i = _LB_ITERS ; --i ; ) {
		jt = Ap.rawBegin();
		it = A.rawBegin();
		for ( ; it != A.rawEnd(); ++it , ++jt) {
			F.init(*jt,*it);
		}
	}
	tim.stop();
	// A.write(std::cout << "--------A---------" << std::endl);
	// Ap.write(std::cout << "--------Ap--------" << std::endl);

	if (!equal(Ap,Bp))
		std::cout << "faux" << std::endl;
	std::cout << "Init  : " << tim << std::endl;
	clear(Ap);
	///////////////////////////////
	tim.clear();
	tim.start();
	for (size_t i = _LB_ITERS ; --i ; ) {
		jt = Ap.rawBegin();
		it = A.rawBegin();
		for ( ; it != A.rawEnd(); ++it , ++jt) {
			*jt = modit(*it,p);
		}
	}
	tim.stop();
	// A.write(std::cout << "--------A---------" << std::endl);
	// Ap.write(std::cout << "--------Ap--------" << std::endl);

	if (!equal(Ap,Bp))
		std::cout << "faux" << std::endl;
	std::cout << "%     : " << tim << std::endl;
	clear(Ap);
	///////////////////////////////
	tim.clear();
	tim.start();
	for (size_t i = _LB_ITERS ; --i ; ) {
		jt = Ap.rawBegin();
		it = A.rawBegin();
		for ( ; it != A.rawEnd(); ++it , ++jt) {
			Integer Jit=*it ; // sinon segfault
			Integer::modin(Jit,(long)p);
			*jt=Jit;
		}
	}
	tim.stop();
	// A.write(std::cout << "--------A---------" << std::endl);
	// Ap.write(std::cout << "--------Ap--------" << std::endl);


	if (!equal(Ap,Bp))
		std::cout << "faux" << std::endl;
	std::cout << "mod   : " << tim << std::endl;
	clear(Ap);

	return 0 ;
}
